In [ ]:
fig = plt.figure()
plt.plot(np.log(np.abs(z_left)),np.log(rho_left), "r-", label = "$z \\in [-L/2,0]$")
plt.plot(np.log(np.abs(z_right)),np.log(rho_right), 'b-', label = "$z \\in [0,L/2]$")
plt.legend()
plt.xlabel("$log|z|$")
plt.ylabel("$log|\\rho|$")
plt.show()

plt.plot(np.log(-z_left[::-1]),np.log(rho_avgd))
plt.legend()
plt.xlabel("$log|z|$")
plt.ylabel("$log|\\rho|$")
plt.show()



#ADDITIONAL:
#Calculate potential 
phi = GF.fourier_potentialV2(rho,L)
#Calculate Acceleration Field on Mesh:
a_grid = NB.acceleration(phi,L) 
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/home/boris/Documents/Research/Coding/1D_Codes/Non-Dim/Analysis/CurveFitting.ipynb Cell 7' in <cell line: 2>()
      <a href='vscode-notebook-cell:/home/boris/Documents/Research/Coding/1D_Codes/Non-Dim/Analysis/CurveFitting.ipynb#ch0000006?line=0'>1</a> fig = plt.figure()
----> <a href='vscode-notebook-cell:/home/boris/Documents/Research/Coding/1D_Codes/Non-Dim/Analysis/CurveFitting.ipynb#ch0000006?line=1'>2</a> plt.plot(np.log(np.abs(z_left)),np.log(rho_left), "r-", label = "$z \\in [-L/2,0]$")
      <a href='vscode-notebook-cell:/home/boris/Documents/Research/Coding/1D_Codes/Non-Dim/Analysis/CurveFitting.ipynb#ch0000006?line=2'>3</a> plt.plot(np.log(np.abs(z_right)),np.log(rho_right), 'b-', label = "$z \\in [0,L/2]$")
      <a href='vscode-notebook-cell:/home/boris/Documents/Research/Coding/1D_Codes/Non-Dim/Analysis/CurveFitting.ipynb#ch0000006?line=3'>4</a> plt.legend()

NameError: name 'z_left' is not defined
<Figure size 432x288 with 0 Axes>
In [ ]:
#hist, xedges, yedges = np.histogram2d(stars_x,stars_v,bins = [50,50])
#plt.imshow(hist,extent = (x_min,x_max,v_min,v_max), cmap = cm.hot, aspect = (x_max-x_min)/(v_max-v_min))
#plt.show()

#Re-center the system
z = np.linspace(-L/2,L/2,N)
max_index = 620
print(z[max_index])
stars_x_new = stars_x - z[max_index] #centroid_z

nbins = 100
hist, xedges, yedges, image = plt.hist2d(stars_x_new, stars_v, 
                            bins = [nbins,nbins], 
                            range = [[-L/2, L/2],[np.min(stars_v),np.max(stars_v)]],
                            #cmax = 250,
                            cmap = cm.hot)
plt.colorbar()
plt.show()

# print(hist)
# for i in range(len(hist)):
#     for j in range(len(hist)):
#         hist[i,j] = int(hist[i,j])
# print(hist)

def Abel(x_range: tuple, v_range: tuple,f):
    # f is a 2d-array. (A 2d histogram)

    x_min,x_max = x_range
    #v_min,v_max = v_range
    
    n_rows,n_cols = np.shape(f)

    #v_s = np.linspace(v_min,v_max,n_rows)
    x_s = np.linspace(x_min,x_max,n_cols)

    #dv = v_s[1]-v_s[0]
    dx = x_s[1]-x_s[0]

    holder = np.array([])
    for i in range(n_rows):
        #v = v_s[i]
        sum = 0
        for j in range(n_cols):
            #x = dx*j #x_s[j]. Note we only add up from x=0 to x->inf
            sum+= f[i][j]*dx
            
            # r = np.sqrt(v**2+x**2)
            # dr = np.sqrt(dv**2+dx**2)
            #
            # term1 = f[i][j]*r
            # term2 = np.sqrt(r**2-v**2)
            # sum+=dr*term1/term2
        holder = np.append(holder,sum)#2*sum)
    return holder

np.savetxt("hist.csv",hist,fmt ='%i',delimiter = ",")
hist = np.loadtxt("hist.csv",dtype = int,delimiter = ',')
v_min,v_max = np.min(stars_v),np.max(stars_v)
abel_transform = Abel((-L/2,L/2),(v_min,v_max),hist)

#print(abel_transform)
v_array = np.linspace(v_min,v_max,len(abel_transform))
plt.plot(v_array,abel_transform,'b--',marker = '.')
plt.show() 
0.24124124124124124
In [ ]:
import Analysis 

#[ ..., [r,m,Num_bosons,sigma,Num_stars],...]
Args = [
    [0.5,1.0,0,1,10000],
    [0.5,1.0,10000,1,10000],
    [1,0.5,20000,1,10000],
    [5,0.1,100000,1,10000],
    [10,0.05,200000,1,10000],
    [50,0.01,1000000,1,10000],
    [0.5,1.0,10000,1,0]
]

for args in Args:
    print("----------New Analysis--------")
    print(
        f"r = {args[0]}",
        f"mu = {args[1]}",
        f"Num_bosons = {args[2]}",
        f"sigma = {args[3]}",
        f"Num_stars = {args[4]}"
    )
    Analysis.analysis(*args)
/home/boris/Documents/Research/Coding/1D_Codes/Non-Dim/Analysis
----------New Analysis--------
r = 0.5 mu = 1.0 Num_bosons = 0 sigma = 1 Num_stars = 10000
/home/boris/Documents/Research/Coding/OneD/WaveNonDim.py:128: RuntimeWarning: invalid value encountered in true_divide
  F_s = F_s/Norm_const
0.002002002002002068
chi^2 = 1069130.1514376372
chi^2 = 331220.8341718659
----------New Analysis--------
r = 0.5 mu = 1.0 Num_bosons = 10000 sigma = 1 Num_stars = 10000
0.002002002002002068
chi^2 = 1923005.6906601514
chi^2 = 1931769.6458091086
chi^2 = 500753.0354785936
----------New Analysis--------
r = 1 mu = 0.5 Num_bosons = 20000 sigma = 1 Num_stars = 10000
0.002002002002002068
chi^2 = 2025205.074171987
chi^2 = 1022913.9619586606
chi^2 = 617435.4881954561
----------New Analysis--------
r = 5 mu = 0.1 Num_bosons = 100000 sigma = 1 Num_stars = 10000
0.002002002002002068
chi^2 = 916046.8264148047
chi^2 = 851932.3919991729
chi^2 = 328193.44588331267
chi^2 = 330860.4024658548
----------New Analysis--------
r = 10 mu = 0.05 Num_bosons = 200000 sigma = 1 Num_stars = 10000
0.0
chi^2 = 1738057.7770754427
chi^2 = 1716889.4310029696
chi^2 = 284629.91746737063
chi^2 = 291832.4249495785
----------New Analysis--------
r = 50 mu = 0.01 Num_bosons = 1000000 sigma = 1 Num_stars = 10000
0.002002002002002068
chi^2 = 309143.8474700038
chi^2 = 321876.3329614536
chi^2 = 494952.2346602648
/home/boris/Documents/Research/Coding/1D_Codes/Non-Dim/Analysis/Analysis.py:437: RuntimeWarning: invalid value encountered in power
  og = a0/((z**a2) * (z+a1)**a3)
chi^2 = 367295.279247727
----------New Analysis--------
r = 0.5 mu = 1.0 Num_bosons = 10000 sigma = 1 Num_stars = 0
In [ ]:
import numpy as np
import matplotlib.pyplot as plt

My_Package_PATH = "/home/boris/Documents/Research/Coding"
import sys
sys.path.insert(1, My_Package_PATH)
import OneD.NBody as NB

z = np.linspace(-1,1)

x = 0
v = 1
star = NB.star(0,1,x,v)

dt = 0.1
t = 0
i = 0
while t < 2:
    plt.plot(star.x,0,'ro')
    plt.xlim(-1,1)
    plt.show()

    star.x -= v*dt
    star.reposition(2)
    t += dt
div = -1.0, remainder = 0.9
new z = -0.09999999999999998
 
div = -1.0, remainder = 0.8
new z = -0.19999999999999996
 
div = -1.0, remainder = 0.7000000000000001
new z = -0.29999999999999993
 
div = -1.0, remainder = 0.6000000000000001
new z = -0.3999999999999999
 
div = -1.0, remainder = 0.5000000000000001
new z = -0.4999999999999999
 
div = -1.0, remainder = 0.40000000000000013
new z = -0.5999999999999999
 
div = -1.0, remainder = 0.30000000000000016
new z = -0.6999999999999998
 
div = -1.0, remainder = 0.20000000000000018
new z = -0.7999999999999998
 
div = -1.0, remainder = 0.1000000000000002
new z = -0.8999999999999998
 
div = -1.0, remainder = 2.220446049250313e-16
new z = -0.9999999999999998
 
div = -2.0, remainder = 0.9000000000000001
new z = 0.9000000000000001
 
div = 0.0, remainder = 0.8000000000000002
new z = 0.8000000000000002
 
div = 0.0, remainder = 0.7000000000000002
new z = 0.7000000000000002
 
div = 0.0, remainder = 0.6000000000000002
new z = 0.6000000000000002
 
div = 0.0, remainder = 0.5000000000000002
new z = 0.5000000000000002
 
div = 0.0, remainder = 0.40000000000000024
new z = 0.40000000000000024
 
div = 0.0, remainder = 0.30000000000000027
new z = 0.30000000000000027
 
div = 0.0, remainder = 0.20000000000000026
new z = 0.20000000000000026
 
div = 0.0, remainder = 0.10000000000000026
new z = 0.10000000000000026
 
div = 0.0, remainder = 2.498001805406602e-16
new z = 2.498001805406602e-16
 
In [ ]: